# PACKAGES
## Tidy code
library(dplyr)
library(magrittr)
library(tibble)
library(purrr)

## Plotting
library(ggplot2)
library(patchwork)
library(reshape2)
library(ggforce)
library(ggpubr)
library(plotly)
library(graphics)
library(circlize)

## Stats
library(MASS)
library(car)
library(FactoMineR)
library(factoextra)
library(ca)
library(vegan)
library(cluster)
library(heplots)
library(rstatix)
library(MVN)

# custom functions
profile_dist_matrix <- function(data){
    mat <- as.matrix((data))
    n <- ncol(mat)
    profiles <- prop.table(x = mat, margin = 2)
    average.profile <- rowSums(mat)/sum(mat)
    dist.mat<- matrix(NA, n, n)
    diag(dist.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            d2 <- sum(((profiles[, i] - profiles[, j])^2) / average.profile)
            dist.mat[i, j] <- dist.mat[j, i] <- d2
        }
    }
  colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
  dist.mat
}

Información

El código y los archivos orginales de esta PEC, así como los informes en formato HTML pueden consultarse en el siguiente repositorio de GitHub. El repositorio es privado, así que se tendrá que requerir el acceso para poder ver el código.

Ejercicios

Ejercicio 1.

En el trabajo de Hunt et al. se estudió la capacidad reproductiva de cinco especies de aves marinas en dos colonias en el sureste del mar de Bering. Además, el apéndice de este estudio resume las colonias y los tamaños de las poblaciones de otros trabajos. El archivo seabirds.csv recoge los datos (número de pájaros) de 23 especies en 9 colonias en el área del norte polar y subpolar. El principal interés de este ejercicio es representar las colonias de diversas formas y estudiar posibles conglomerados.

seabirds <- read.csv("seabirds.csv", stringsAsFactors = F)
seabirds2 <- seabirds %>% column_to_rownames("Specie")
margin_col <- rowSums(seabirds2)
margin_row <- colSums(seabirds2)
n <- sum(seabirds2)

A)

Calcular las frecuencias relativas, las frecuencias relativas marginales y la matriz de perfiles. El resultado debería ser la tabla 12.6 del libro de Krebs y que reproducimos al final de este documento.

Para calcular las frecuencias relativas, solo tenemos que dividir cada valor (frecuencias absolutas) por el total de observaciones (4765889) o usando la función prop.table(x=seabirds) y las frecuencias relativas marginales se consiguen con la suma de las frecuencias relativas de cada columna o de cada fila. Sin embargo, con esto no nos sale igual que la tabla 12.6 del libro de Krebs.

seabirds_freq_rel2 <- prop.table(x = seabirds2 %>% as.matrix())
addmargins(seabirds_freq_rel2) %>% round(4) %>% 
  knitr::kable(caption = "Table 1.1. Relative and marginal frequencies from data in 'seabirds.csv'. Relative frequencies are calculated using the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name 'sum'. This table is not coincident with table 12.6 from Krebs book.", align = "c")
Table 1.1. Relative and marginal frequencies from data in ‘seabirds.csv’. Relative frequencies are calculated using the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name ‘sum’. This table is not coincident with table 12.6 from Krebs book.
CH PLI CI NS CL CT SI SPI SGI Sum
Northern.fulmar 0.0000 0.0260 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0147 0.0409
Glaucous-winged.gull 0.0000 0.0001 0.0000 0.0001 0.0000 0.0001 0.0000 0.0000 0.0000 0.0003
Black-legged.kittiwake 0.0084 0.0122 0.0126 0.0016 0.0052 0.0056 0.0010 0.0065 0.0151 0.0682
Red-legged.kittiwake 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0005 0.0462 0.0466
Thick-billed.murre 0.0588 0.0361 0.0671 0.0001 0.0063 0.0490 0.0000 0.0231 0.3147 0.5552
Common.murre 0.0000 0.0000 0.0000 0.0088 0.0147 0.0327 0.0011 0.0082 0.0399 0.1053
Black.guillemot 0.0000 0.0017 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0018
Pigeon.guillemot 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Horned.puffin 0.0000 0.0000 0.0000 0.0007 0.0003 0.0003 0.0000 0.0009 0.0059 0.0081
Tufted.puffin 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0013 0.0015
Atlantic.puffin 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0034 0.0000 0.0000 0.0034
Pelagic.cormorant 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001
Red-faced.cormorant 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0005 0.0010 0.0016
Shag 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Parakeet.auklet 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0071 0.0315 0.0386
Crested.auklet 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0013 0.0059 0.0071
Least.auklet 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0048 0.0525 0.0573
Razorbill 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.0000 0.0000 0.0009
Manx.shearwater 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0546 0.0000 0.0000 0.0546
Storm.petrel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0027 0.0000 0.0000 0.0027
Herring.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0016 0.0000 0.0000 0.0016
Great.black-backed.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0001
Lesser.black.backed.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0042 0.0000 0.0000 0.0042
Sum 0.0672 0.0760 0.0798 0.0113 0.0266 0.0876 0.0696 0.0533 0.5285 1.0000

Para que coincida con la tabla 12.6, la tabla de frecuencias relativas debe hacerse respecto al total de observaciones en cada lugar, por lo que el valor de la frecuencia absoluta dividido por la suma de todas las observaciones en la columna correspondiente, algo que también pude conseguirse con prop.table(x = seabirds, margin = 2). Lo que es la matriz de perfiles de las columnas o la frecuencia relativa condicionada a las columnas (frec. absoluta de una celda dividida por la suma de las frec. absolutas de la columna correspondiente).

seabirds_perfiles_c <- prop.table(x = seabirds2 %>% as.matrix(), margin = 2)
addmargins(seabirds_perfiles_c) %>% round(4) %>% 
  knitr::kable(caption = "Table 1.2. Relative and marginal frequencies conditioned to columns from data in 'seabirds.csv'. Relative frequencies are calculated taking all the observations in each place (column), not the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name 'sum'.", align = "c")
Table 1.2. Relative and marginal frequencies conditioned to columns from data in ‘seabirds.csv’. Relative frequencies are calculated taking all the observations in each place (column), not the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name ‘sum’.
CH PLI CI NS CL CT SI SPI SGI Sum
Northern.fulmar 0.0000 0.3422 0.0000 0.0000 0.0000 0.0000 0.0007 0.0028 0.0278 0.3734
Glaucous-winged.gull 0.0005 0.0011 0.0004 0.0051 0.0004 0.0007 0.0000 0.0000 0.0000 0.0082
Black-legged.kittiwake 0.1249 0.1600 0.1577 0.1402 0.1972 0.0634 0.0151 0.1221 0.0286 1.0094
Red-legged.kittiwake 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0087 0.0873 0.0960
Thick-billed.murre 0.8740 0.4746 0.8413 0.0074 0.2367 0.5593 0.0000 0.4334 0.5955 4.0222
Common.murre 0.0000 0.0000 0.0000 0.7765 0.5522 0.3728 0.0160 0.1537 0.0754 1.9466
Black.guillemot 0.0006 0.0221 0.0005 0.0000 0.0013 0.0000 0.0000 0.0000 0.0000 0.0246
Pigeon.guillemot 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Horned.puffin 0.0000 0.0000 0.0000 0.0592 0.0114 0.0036 0.0000 0.0173 0.0111 0.1026
Tufted.puffin 0.0000 0.0000 0.0000 0.0008 0.0002 0.0000 0.0000 0.0039 0.0024 0.0072
Atlantic.puffin 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0482 0.0000 0.0000 0.0482
Pelagic.cormorant 0.0000 0.0000 0.0000 0.0096 0.0006 0.0001 0.0001 0.0000 0.0000 0.0104
Red-faced.cormorant 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0099 0.0020 0.0118
Shag 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0001
Parakeet.auklet 0.0000 0.0000 0.0000 0.0012 0.0000 0.0000 0.0000 0.1340 0.0595 0.1947
Crested.auklet 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0236 0.0111 0.0348
Least.auklet 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0906 0.0992 0.1899
Razorbill 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0130 0.0000 0.0000 0.0130
Manx.shearwater 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.7838 0.0000 0.0000 0.7838
Storm.petrel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0389 0.0000 0.0000 0.0389
Herring.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0229 0.0000 0.0000 0.0229
Great.black-backed.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.0000 0.0000 0.0009
Lesser.black.backed.gull 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0603 0.0000 0.0000 0.0603
Sum 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 9.0000

B)

Calcular la matriz de distancias ji-cuadrado entre los perfiles de las columnas y su inercia total.

Los perfiles de las columnas estan calculados en el ejercicio enterior, así que usaremos esos datos para calcular las distancias ji-cuadrado y su inercia total.

Distancia ji-cuadrado: con la distancia ji-cuadrado entre los perfiles, podemos observar qué columnas se parecen más unas a otras: a una distancia más pequeña, más similitud entre columnas (en este caso, colonias de varias aves).

La distancia ji-cuadrado entre columnas se puede definir como

\[ d_{ij}^{(cols)} = \sum_{k=1}^{r} \frac{1}{p_k.}(p^c_{ki} - p^c_{kj})^2 \]

donde:

  • \(d_{ij}^{(cols)}\) es la distancia entre dos columnas \(i\) y \(j\),
  • \(r\) es el número de filas,
  • \(p_k. = \frac{n_k.}{n.}\) con \(n_k.\) siendo la frecuencia marginal de la fila \(k\) y \(n.\) el total de observaciones de la muestra (también dicho el perfil de columnas medio),
  • \(p^c_{ki}\) es la proporción en la fila \(k\) y la columna \(i\) y
  • \(p^c_{kj}\) es la proporción en la fila \(k\) y la columna \(j\)

En nuestro caso, las proporciones las tenemos en la tabla calculada en el ejercicio anterior (ejercicio 1a) y que coinciden con la tabla 12.6 del libro de Krebs.

Para calcular dicha matriz, usaremos la función profile_dist_matrix(), definida a l’inicio de este informe, sacada de STHDA y modificada por mi. profile_dist_matrix() acepta como input la matriz de perfiles (que devemos de transformar con t() para que coja las columnas) y el perfil de columnas medio (\(p_k.\) en la función de la distancia entre columnas)

dist_cols_matrix <- profile_dist_matrix((seabirds2))
dist_cols_matrix %>% knitr::kable(caption = "Table 1.3. Chi-squared distance matrix between column profiles of the 'seabirds.csv' data.")
Table 1.3. Chi-squared distance matrix between column profiles of the ‘seabirds.csv’ data.
CH PLI CI NS CL CT SI SPI SGI
CH 0.0000000 3.425238 0.0178094 8.251115 3.7234212 1.5564260 15.46510 1.3725624 0.8159060
PLI 3.4252376 0.000000 3.3647355 10.407262 6.1401304 4.6055964 17.75168 4.1374664 3.4917128
CI 0.0178094 3.364736 0.0000000 8.154352 3.5963874 1.5964295 15.48546 1.3407506 0.8934179
NS 8.2511147 10.407262 8.1543518 0.000000 1.5595341 3.2777270 20.80417 5.7319920 6.9746456
CL 3.7234212 6.140130 3.5963874 1.559534 0.0000000 0.7664893 17.24585 2.4287149 3.2811626
CT 1.5564260 4.605596 1.5964295 3.277727 0.7664893 0.0000000 15.71938 1.3213094 1.3394868
SI 15.4650970 17.751675 15.4854630 20.804173 17.2458462 15.7193804 0.00000 15.3926081 15.0680100
SPI 1.3725624 4.137466 1.3407506 5.731992 2.4287149 1.3213094 15.39261 0.0000000 0.5942528
SGI 0.8159060 3.491713 0.8934179 6.974646 3.2811626 1.3394868 15.06801 0.5942528 0.0000000

Como és de esperar, en la diagonal observamos distancias de 0, ya que comparamos cada columna consigo misma.

Por otro lado, podemos ver como hay columnas con una distancia muy pequeña como CH-CI (0.0178) o CL-CT (0.7665), mientras que otras tienen una distancia mayor y, por lo tanto, observamos como són distintas (CH-SI = 15.4651, NS-SI = 20.8042).

Inercia total: la inercia total es la media ponderada de los cuadrados de las distancias \(\chi^2\) entre los perfiles fila/columna y el perfil media [^(2)^](https://www.fbbva.es/wp-content/uploads/2017/05/dat/greenacre_cap04.pdf). También puede encontrarse dividiendo el estadístico \(\chi^2\) por el total de observaciones. Cabe destacar que la inercia total es la misma tanto si se hace con los perfiles de fila como los de columna, ya que estamos actuando sobre la totalidad de los datos.

as.numeric(chisq.test(seabirds2)$stat)/n
## [1] 1.565692

C)

Con la matriz de distancias ji-cuadrado entre los perfiles realizar un escalado multidimensional. Dibujar las coordenadas principales para las columnas.

Para el escalado multidimensional usaremos la función cmdscale() con 8 coordenadas principales (el máximo en este caso) y indicando que queremos que nos devuelva los valores propios (eigenvalues).

c <- cmdscale(d = dist_cols_matrix, k=8, eig = T)
c
## $points
##            [,1]        [,2]        [,3]        [,4]         [,5]
## CH    0.5043776 -1.75638463 -0.72114533  0.02370077  0.120329377
## PLI   1.9436165 -4.95136344  0.48741897  0.10467367 -0.010646143
## CI    0.5461064 -1.70615994 -0.70169181 -0.09456648 -0.110423324
## NS    6.0813689  4.66506492  0.15933028  0.03698951  0.001837342
## CL    2.8046648  1.78390484 -0.31619954 -0.11640599 -0.016327900
## CT    1.4177905  1.07090468  0.05751178  0.48401878 -0.004345404
## SI  -14.5809325  1.59022765  0.06534089  0.02022302 -0.002890291
## SPI   0.8713107  0.05057765  0.41238536 -0.57607534  0.021088565
## SGI   0.4116971 -0.74677172  0.55704940  0.11744207  0.001377778
## 
## $eig
## [1]  2.647219e+02  5.969290e+01  1.863310e+00  6.157183e-01  2.752965e-02
## [6] -3.197442e-14 -6.373612e-01 -1.598085e+00 -1.998594e+01
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 0.9363544 1.0000000

Una vez realizado el escalado multidimensional, debemos ver cuantas dimensiones (o coordenadas principales) explican la mayor parte de los datos originales.

Para ello usaremos los dos métodos siguientes, que se basan en la suma acumulada de los valores propios dividida por la suma total de estos valores propios. Como en nuestro caso hay algunos valores negativos, es necesario usar el valor absoluto (método 1) o elevar al cuadrado (método 2).

\[ P_{m}^{(1)} = \frac{\sum^{m}_{i=1} |\lambda_{i}| }{\sum^{n}_{i=1} |\lambda_{i}| } \]

cumsum(abs(c$eig)) / sum(abs(c$eig))
## [1] 0.7582053 0.9291753 0.9345121 0.9362756 0.9363544 0.9363544 0.9381799
## [8] 0.9427571 1.0000000

\[ P_{m}^{(1)} = \frac{\sum^{m}_{i=1} \lambda_{i}^2 }{\sum^{n}_{i=1} \lambda_{i}^2 } \]

cumsum(c$eig^2) / sum(c$eig^2)
## [1] 0.9463924 0.9945136 0.9945605 0.9945656 0.9945657 0.9945657 0.9945711
## [8] 0.9946056 1.0000000

Ambos métodos nos sugieren utilizar dos dimensiones para la representación de los datos, pues después de la segunda componente principal, las otras no parecen tener un gran papel en la explicación de los datos

c$points %>% as.data.frame() %>% mutate(names = rownames(c$points)) %>%
ggplot(aes(V1, V2)) +
  geom_point() +
  geom_text(aes(label = names), 
            position = position_dodge(width = 1),
            vjust = -1, hjust = 1, size = 4) +

  ggtitle("MDS of Chi-squared distances",
          "of column profiles") +
  xlab("1st principal coordinate") + ylab("2nd principal coordinate") +
  
  theme_bw(base_size = 15)
<b>Figure 1.1. MDS of chi-squared distances of column profiles. 2D representation (dimensións 1 and 2).</b>

Figure 1.1. MDS of chi-squared distances of column profiles. 2D representation (dimensións 1 and 2).

Como se puede ver, la mayoría de columnas se encuentran en un punto similar dentro de la primera coordinada principal (salvo la columna SI), que es la que explica la mayor parte de la variación, mientras que en la segunda coordinada varian mucho.

Vemos como con dos dimensiones se explica mucho la distancia/dissimilaridad entre columnas, ya que cuanto más alejadas están unas de otras en el gráfico de coordenadas principales, más distancia hay entre sus perfiles. Por ejemplo, SI vemos como tiene un perfil que tiene una distáncia grande con todas las otras, PLI y NS también se encuentran lejos entre sí, mientras que CI y CH están muy cerca.

A continuación, aunque no es extrictamente necesario, vamos ha hacer un gráfico 3D con las primeras tres coordinadas principales.

plot3d <- cmdscale(dist_cols_matrix, k=8) %>% as.data.frame() %>% rownames_to_column("Name")
plotly::plot_ly(x = plot3d$V1, y = plot3d$V2, z = plot3d$V3, name = plot3d$Name,
                type = "scatter3d")

Figura 1.2. MDS of chi-squared distances of column profiles. 3D representation (dimensións 1, 2 and 3).

Podemos ver como en el eje Z (tercera coordinada principal), las columnas no parecen tener mucha variabilidad, dado que el rango de Z es de 1, mientras que los otros dos son de 20 y 8 para X e Y respectivamente.

Aún así, vemos como esta nueva coordenada principal explica un poco más la distancia entre los perfiles de las columnas. Por ejemplo, las columnas CI y CH, que con dos dimensiones se veían muy cerca de SGI, ahora se encuentran en los extremos opuestos del eje Z.

Sin embargo, esta representación puede ser confusa, ya que debido a las escalas del gráficovemos como, por ejemplo, CH se encuentra muy distanciada de SGI y parece estar más cerca de CL, cuando las distancias ji-cuadrado entre sus perfiles indican lo contrario.

D)

Realizar un análisis de correspondencias y calcular las inercias principales (en %) y la inercia total con los valores propios. Dibujar una representación simétrica del CA. A pesar de la confusión de nombres, ¿cuales son las especies que caracterizan a la colonia SI (Skomer Island, Irish Sea)?

El análisis de correspondencias vamos a hacerlo usando la función CA() del paquete FactoMineR. Esta función nos devuelve los valores propios (eigenvalues) y los porcentajes de variancia (inercia) de cada dimensión ademas de las coordenadas de las filas y las columnas en cada dimensión y sus inercias, contribuciones, etc. Asimismo, nos permite hacer directamente una representación de los datos (indicando graph=T).

Como el porcentaje de inercia ya nos lo da la función, no hace falta calcularlo, aunque sencillamente sería calcular la suma de los valores propios (inercia total) y efectuar el porcentaje de cada valor propio respecto a ese valor.

\[ Inercia \space total = {\sum_{i=1}^{n}\lambda_i}*100 \] Siendo \(\lambda_i\) los valores propios y \(n\) el total de dimensiones del anàlisis de correspondencias.

\[ Inercia \space principal \space _{(de \space la \space dimension \space j)} = \frac{\lambda_j}{\sum_{i=1}^n\lambda_i}*100 \]

ca = CA(X = seabirds2, ncp = 5, graph = F)

Los valores propios, así como la inercia principal (en porcentaje) y el porcentaje cumulativo de variancia/inercia se encuentran dentro la el resultado de la funcón CA(), concretamente en la variable $eig.

ca$eig %>% knitr::kable(caption = "Table 1.4. Eigenvalues, percentages of variance/inertia and cumulative percentage of variance/inertia for each dimention.")
Table 1.4. Eigenvalues, percentages of variance/inertia and cumulative percentage of variance/inertia for each dimention.
eigenvalue percentage of variance cumulative percentage of variance
dim 1 0.9662498 61.7139287 61.71393
dim 2 0.2566386 16.3913870 78.10532
dim 3 0.2059249 13.1523299 91.25765
dim 4 0.0986173 6.2986387 97.55628
dim 5 0.0261770 1.6719151 99.22820
dim 6 0.0083129 0.5309395 99.75914
dim 7 0.0037650 0.2404718 99.99961
dim 8 0.0000061 0.0003892 100.00000

En los datos anteriores podemos ver los distintos porcentajes de inercia de cada dimensión y se ve claramente que la dimensión 1 nos explica la mayor parte de la variación (\(\pm\) 61.7%), mientras que las dimensiones 2 y 3 (\(\pm\) 16.4 y \(\pm\) 13.2, respectivamente) nos explican un porcentaje relativamente mayor que las siguientes ddimensiones, que nos explican muy poco en total.

Para calcular la inercia total tenemos que sumar todos los valores propios, tal y como hemos mencionado anteriormente.

paste("La inercia total es:", sum(ca$eig[,1]))
## [1] "La inercia total es: 1.56569159050875"

Finalmente, para representar los datos podemos usar las coordenadas de las filas y las columnas en las dos primeras dimensiones y hacer un plot() de ellas o simplemente indicando graph=T dentro de la función CA() o usando ggplot2.

ggplot() + 
  geom_point(ca$row$coord %>% as.data.frame(), mapping = aes(x = `Dim 1`, y = `Dim 2`), color="Blue") +
  ggrepel::geom_text_repel(ca$row$coord %>% as.data.frame() %>% rownames_to_column("Name"), 
            mapping = aes(x = `Dim 1`, y = `Dim 2`, label = Name), color = "Darkblue", seed = 3) +
  
  geom_point(ca$col$coord %>% as.data.frame(), mapping = aes(x = `Dim 1`, y = `Dim 2`), color="Red") +
  ggrepel::geom_label_repel(ca$col$coord %>% as.data.frame() %>% rownames_to_column("Name"), 
            mapping = aes(x = `Dim 1`, y = `Dim 2`, label = Name), color = "Darkred", seed = 3) +
  
  theme_bw(base_size = 14) +
  
  ggtitle("Representation of correspondence analysis",
          "Dimension 1 and dimension 2") +
  xlab("Dimension 1 (61.7% of variance)") +
  ylab("Dimension 2 (16.4% of variance)") 
<b>Figure 1.3. CA representation (dimensions 1 and 2) from 'seabirds.csv', with rows (species) represented in blue and columns (colonies) represented in red.<b>

Figure 1.3. CA representation (dimensions 1 and 2) from ‘seabirds.csv’, with rows (species) represented in blue and columns (colonies) represented in red.

Como se puede ver en la figura 1.3, hay muchas especies características de la colonia SI, entre las cuales se encuentran:

  • Razorbill
  • Herring gull
  • Lesser black backed gull
  • Great black backed gull
  • Manx shearwater
  • Shag
  • Storm petrel
  • Atlantic puffin

E)

E.1.

Dada la gran cantidad de ceros en la tabla 12.6, en el libro de Krebs se sugiere la utilización de la distancia de Canberra entre las columnas de la tabla 12.6. La distancia de Canberra no tiene una única definición y, además, ha cambiado a lo largo de la historia. Una posible definición entre dos vectores \(\textbf{p} = (p_1, p_2, . . . , p_k)'\) y \(\textbf{q} = (q_1, q_2, . . . , q_k)'\) de la misma longitud es

\[ d_C(\textbf{p},\textbf{q}) = \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]

Cuando el denominador es cero, el cociente es NaN, y el sumando se elimina.

Comprobar que esta definición no sirve para calcular la matriz de similaridades de la tabla 12.7 del libro de Krebs y que se reproduce al final de este documento.

mat <- as.matrix(seabirds2)
n <- ncol(mat)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0

for(i in 1:(n-1)){
  for(j in (i+1):n){
    dc <- sum(abs(profiles[, i] - profiles[, j]) /
                (abs(profiles[, i]) + abs(profiles[, j])), na.rm = T)
    dist.mat[i, j] <- dist.mat[j, i] <- dc
  }
}

colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat %>% knitr::kable(caption = "Table 1.5. Canberra distance between column profiles using the 1st definition (mentioned above): d(p,q) = sum(abs(p-q) / (abs(p)+abs(q))), where p and q are 2 column profiles.")
Table 1.5. Canberra distance between column profiles using the 1st definition (mentioned above): d(p,q) = sum(abs(p-q) / (abs(p)+abs(q))), where p and q are 2 column profiles.
CH PLI CI NS CL CT SI SPI SGI
CH 0.0000000 2.768917 0.3066923 7.872672 5.248740 5.726277 14.78458 11.34803 11.81702
PLI 2.7689165 0.000000 2.7128716 8.678678 6.797518 6.724761 14.82392 11.16369 11.65965
CI 0.3066923 2.712872 0.0000000 7.897405 5.108834 5.883052 14.82557 11.44721 11.86432
NS 7.8726718 8.678678 7.8974045 0.000000 6.345527 8.318004 16.74805 10.90990 11.61912
CL 5.2487396 6.797518 5.1088341 6.345527 0.000000 5.609152 16.56029 11.22128 11.82774
CT 5.7262773 6.724761 5.8830525 8.318004 5.609152 0.000000 15.66589 12.51723 12.58676
SI 14.7845812 14.823919 14.8255729 16.748052 16.560290 15.665893 0.00000 19.19796 18.91136
SPI 11.3480265 11.163693 11.4472135 10.909901 11.221283 12.517229 19.19796 0.00000 4.67864
SGI 11.8170207 11.659654 11.8643193 11.619124 11.827739 12.586763 18.91136 4.67864 0.00000

La tabla 12.7 del libro de Krebs muestra las similaridades entre columnas que se han obtenido calculando \(1 - Distancia\space de \space Canberra\). Como vemos en la tabla 1.5, esta muestra distancias mayores que 1, por lo que restarlas a 1 resultaría en similaridades negativas, cosa que no se muestra en la tabla 12.7 del libro de Krebs.

E.2.

Una modificación de la distancia anterior es considerar la distancia

\[ d_C(\textbf{p},\textbf{q}) = \frac{1}{k} \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]

Igual que antes, cuando el denominador es cero, el sumando se elimina.

Comprobar que con esta definición se puede obtener la tabla 12.7.

mat <- as.matrix(seabirds2)
n <- ncol(mat)
k <- nrow(mat)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0

for(i in 1:(n-1)){
  for(j in (i+1):n){
    dc <- sum(abs(profiles[, i] - profiles[, j]) /
              (abs(profiles[, i]) + abs(profiles[, j])), na.rm = T) / k
      
    dist.mat[i, j] <- dist.mat[j, i] <- dc
  }
}

colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)

dist.mat %>% knitr::kable(caption = "Table 1.6. Canberra distance between column profiles using the 2nd definition (mentioned above): d(p,q) = (1/k)*(sum(abs(p-q) / (abs(p)+abs(q)))), where p and q are 2 column profiles and k is the number of rows in the profiles table.")
Table 1.6. Canberra distance between column profiles using the 2nd definition (mentioned above): d(p,q) = (1/k)*(sum(abs(p-q) / (abs(p)+abs(q)))), where p and q are 2 column profiles and k is the number of rows in the profiles table.
CH PLI CI NS CL CT SI SPI SGI
CH 0.0000000 0.1203877 0.0133344 0.3422901 0.2282061 0.2489686 0.6428079 0.4933925 0.5137835
PLI 0.1203877 0.0000000 0.1179509 0.3773338 0.2955442 0.2923809 0.6445182 0.4853779 0.5069415
CI 0.0133344 0.1179509 0.0000000 0.3433654 0.2221232 0.2557849 0.6445901 0.4977049 0.5158400
NS 0.3422901 0.3773338 0.3433654 0.0000000 0.2758925 0.3616524 0.7281762 0.4743435 0.5051793
CL 0.2282061 0.2955442 0.2221232 0.2758925 0.0000000 0.2438762 0.7200126 0.4878819 0.5142495
CT 0.2489686 0.2923809 0.2557849 0.3616524 0.2438762 0.0000000 0.6811258 0.5442274 0.5472506
SI 0.6428079 0.6445182 0.6445901 0.7281762 0.7200126 0.6811258 0.0000000 0.8346939 0.8222330
SPI 0.4933925 0.4853779 0.4977049 0.4743435 0.4878819 0.5442274 0.8346939 0.0000000 0.2034191
SGI 0.5137835 0.5069415 0.5158400 0.5051793 0.5142495 0.5472506 0.8222330 0.2034191 0.0000000

La tabla 1.6 no muestra ninguna distancia mayor que uno, por lo que es posible que, si restamos estas distancias a 1, el resultado sea la tabla 12.7 del libro de Krebs.

round(1-dist.mat, 2) %>% knitr::kable(caption = "Table 1.7. Table of similarities obtained by resting the Canberra distances in table 1.6 from 1. As it can be seen, this table is the same as the table 12.7 from the Krebs' book.")
Table 1.7. Table of similarities obtained by resting the Canberra distances in table 1.6 from 1. As it can be seen, this table is the same as the table 12.7 from the Krebs’ book.
CH PLI CI NS CL CT SI SPI SGI
CH 1.00 0.88 0.99 0.66 0.77 0.75 0.36 0.51 0.49
PLI 0.88 1.00 0.88 0.62 0.70 0.71 0.36 0.51 0.49
CI 0.99 0.88 1.00 0.66 0.78 0.74 0.36 0.50 0.48
NS 0.66 0.62 0.66 1.00 0.72 0.64 0.27 0.53 0.49
CL 0.77 0.70 0.78 0.72 1.00 0.76 0.28 0.51 0.49
CT 0.75 0.71 0.74 0.64 0.76 1.00 0.32 0.46 0.45
SI 0.36 0.36 0.36 0.27 0.28 0.32 1.00 0.17 0.18
SPI 0.51 0.51 0.50 0.53 0.51 0.46 0.17 1.00 0.80
SGI 0.49 0.49 0.48 0.49 0.49 0.45 0.18 0.80 1.00

Como se puede ver en la tabla 1.7, las similaridades son las mismas que en la tabla 12.7 del libro de Krebs (aunque la tabla 1.7 es simetrica respecto la diagonal y la del libro de krebs solo muestra la parte superior de esta diagonal). Por lo tanto, podemos afirmar que la definición usada en la tabla 1.6 para calcular las distancias de Canberra es la misma que la que usan los autores de la 12.7 del libro de Krebs.

E.3.

Finalmente, se puede comprobar que ésta tampoco es la definición que utiliza R para calcular la distancia de Canberra. Tras una ardua investigación, se comprueba que la definición de R es

\[ d_C(\textbf{p},\textbf{q}) = \frac{k}{k-n_z} \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]

donde nz es el número de denominadores cero.

Comprobar que ésta es la definición de la distancia de Canberra según R.

Según la documentación de la función dist(), el cálculo de la distancia de Canberra usando la función dist(x, method = "canberra") del paquete stats no es la que se menciona arriba, sinó la primera (vease la [documentación dedist()`](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/dist)).

Sin embargo, si usamos la definición mencionada en este ejercicio 1.e.3 y calculamos la distancia de Canberra (tabla 1.8) obtenemos las mismas distancias que usando la función dist() de R (tabla 1.9), por lo que podemos afirmar que R usa la definicion mencionada en este ejercicio 1.e.3.

mat <- as.matrix(seabirds2)
n <- ncol(mat)
k <- nrow(mat)
#nz <- sum(mat==0)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0

for(i in 1:(n-1)){
  for(j in (i+1):n){
    nz <- sum((abs(profiles[, i]) + abs(profiles[, j]))==0)
    dc <- sum(abs(profiles[, i] - profiles[, j]) /
              (abs(profiles[, i]) + abs(profiles[, j])), na.rm = T) * (k/(k-nz))
      
    dist.mat[i, j] <- dist.mat[j, i] <- dc
  }
}

colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat %>% 
  knitr::kable(caption = "Table 1.8. Canberra distance calculated using the column profiles of 'seabirds.csv' dataset and the definition given in the exercise 1.E.3: d(p,q) = (k/(k-nz)*(sum(abs(p - q) / (abs(p) + abs(q))))")

(dist_canb <- dist(t(profiles), method = "canberra", diag = T, upper = T)) %>% as.matrix() %>%
  knitr::kable(caption = "Table 1.9. Canberra distance calculated using the column profiles of 'seabirds.csv' dataset and the 'dist()' function specifying `method='canberra'`")
Table 1.8. Canberra distance calculated using the column profiles of ‘seabirds.csv’ dataset and the definition given in the exercise 1.E.3: d(p,q) = (k/(k-nz)*(sum(abs(p - q) / (abs(p) + abs(q))))
CH PLI CI NS CL CT SI SPI SGI
CH 0.000000 12.73702 1.763481 20.11905 15.09013 16.46305 22.66969 20.077278 20.907037
PLI 12.737016 0.00000 12.479210 19.96096 17.37143 17.18550 22.73001 19.751149 20.628619
CI 1.763481 12.47921 0.000000 20.18226 14.68790 16.91378 22.73255 20.252762 20.990719
NS 20.119050 19.96096 20.182256 0.00000 16.21635 19.13141 22.65913 19.302132 20.556911
CL 15.090126 17.37143 14.687898 16.21635 0.00000 14.33450 22.40510 18.434965 19.431285
CT 16.463047 17.18550 16.913776 19.13141 14.33450 0.00000 21.19503 19.193085 19.299703
SI 22.669691 22.73001 22.732545 22.65913 22.40510 21.19503 0.00000 22.077655 21.748062
SPI 20.077278 19.75115 20.252762 19.30213 18.43497 19.19308 22.07765 0.000000 9.782611
SGI 20.907037 20.62862 20.990719 20.55691 19.43129 19.29970 21.74806 9.782611 0.000000
Table 1.9. Canberra distance calculated using the column profiles of ‘seabirds.csv’ dataset and the ‘dist()’ function specifying method='canberra'
CH PLI CI NS CL CT SI SPI SGI
CH 0.000000 12.73702 1.763481 20.11905 15.09013 16.46305 22.66969 20.077278 20.907037
PLI 12.737016 0.00000 12.479210 19.96096 17.37143 17.18550 22.73001 19.751149 20.628619
CI 1.763481 12.47921 0.000000 20.18226 14.68790 16.91378 22.73255 20.252762 20.990719
NS 20.119050 19.96096 20.182256 0.00000 16.21635 19.13141 22.65913 19.302132 20.556911
CL 15.090126 17.37143 14.687898 16.21635 0.00000 14.33450 22.40510 18.434965 19.431285
CT 16.463047 17.18550 16.913776 19.13141 14.33450 0.00000 21.19503 19.193085 19.299703
SI 22.669691 22.73001 22.732545 22.65913 22.40510 21.19503 0.00000 22.077655 21.748062
SPI 20.077278 19.75115 20.252762 19.30213 18.43497 19.19308 22.07765 0.000000 9.782611
SGI 20.907037 20.62862 20.990719 20.55691 19.43129 19.29970 21.74806 9.782611 0.000000

F)

Realizar un MDS con la distancia de Canberra de R. Comprobar que se trata de una distancia euclídea. Dibujar el mapa. Comparar el resultado con el obtenido con la distancia ji-cuadrado. Utilizar la función procrustes() del paquete vegan.

mds2 <- cmdscale(d = dist_canb, k = 8, eig = T)
mds2$points %>% knitr::kable(caption = "Table 1.10. Coordinates of the MDS from Canberra distances between column profiles shown in table 1.9.")
mds2$eig %>% knitr::kable(caption = "Table 1.11. Eigenvalues of the MDS from Canberra distances between column profiles shown in the table 1.9.")
Table 1.10. Coordinates of the MDS from Canberra distances between column profiles shown in table 1.9.
CH -9.004001 -0.4119965 3.403355 -1.3822228 -2.9394413 2.7220433 0.2941079 0.6508961
PLI -6.548694 -0.5054878 3.377353 -3.1878164 7.5606813 -3.7432170 -0.4103127 0.0303465
CI -9.199724 -0.4306546 3.359156 -1.6881892 -3.2046800 1.7311915 -0.1130372 -0.6748997
NS 2.405435 -0.5506840 -11.880669 -6.1205145 1.2539361 2.3449580 -0.3070384 -0.0131285
CL -2.708016 -1.5589356 -5.787390 3.8178859 -4.6417387 -5.7230153 -0.5299328 0.0591010
CT -1.427421 1.0272123 -2.678522 9.8718390 3.9359064 2.9963500 0.5743877 -0.0481728
SI 7.209240 14.9136546 2.600906 -0.9028339 -1.0232724 -0.6245509 0.2122803 0.0041246
SPI 9.087292 -6.5746806 3.163752 -0.8795986 -0.5112311 -0.4970329 4.8545370 -0.0253470
SGI 10.185890 -5.9084277 4.442059 0.4714505 -0.4301603 0.7932733 -4.5749918 0.0170798
Table 1.11. Eigenvalues of the MDS from Canberra distances between column profiles shown in the table 1.9.
x
462.0537185
304.9524572
252.5978579
166.2239203
116.1771250
72.9144505
45.5147980
0.8870133
0.0000000
mds2$points %>% as.data.frame() %>% mutate(names = rownames(c$points)) %>%
ggplot(aes(V1, V2)) +
  geom_point() +
  geom_text(aes(label = names), 
            position = position_dodge(width = 1),
            vjust = -1, hjust = 1, size = 4) +

  ggtitle("MDS of Canberra distances",
          "of column profiles") +
  xlab("1st principal coordinate") + ylab("2nd principal coordinate") +
  
  theme_bw(base_size = 15)
<b>Figure 1.4. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)<b>

Figure 1.4. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)

plot(procrustes(c$points, mds2$points, scale = T))
<b>Figure 1.5. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)<b>

Figure 1.5. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)

Dos características que tienen las matrices de distancias no Euclídeas es que algunos de los valores propios del escalado multidimensional seran negativos, mientras que algunas de las coordenadas pueden ser numeros complejos. En nuestro caso, ni las coordenadas (tabla 1.10) son numeros complejos, ni los valores propios (tabla 1.11) son negativos. Por este motivo podemos sugerir que las distancias de la tabla 1.9 son Euclídeas.

G)

Realizar un análisis de conglomerados jerárquico con el método deWard2 de la distancia de Canberra según R. Dibujar el dendograma resultante. Dibujar también un heatmap de este análisis con la función heatmap(). Para ello, hay que elegir bien los parámetros distfun= y hclustfun= y una escala de colores. ¿Para qué sirve el heatmap?

hclust_1g <- hcut(x = dist_canb, hc_method = "ward.D2", k = 6)
fviz_dend(hclust_1g, cex = 2 , k = 6, main = "Hierarchical clustering with Canberra distances", sub = "6 clusters",
          ggtheme = theme_bw(base_size = 17), show_labels = T, rect = TRUE)
<b>Figure 1.6.</b>

Figure 1.6.

Con el análisis de conglomerados jerárquico (figura 1.6), observamos como muchas de las afirmaciones que hicimos en análisis anteriores salen reforzadas. Para analizar el dendograma de la figura 1.6, debemos observar a que altura se separan cada conglomerado, hasta llegar a los conglomerados (clusters) donde se encuentran las columnas individuales. Cuanto a más altura se de una separción, más lejos estaran los clusteres resultatnes y, para comparar dos clusters, podemos sumar las alturas que los separan:

  • Primero vemos como hay una separación inicial que agrupa SI, SPI y SGI en un conglomerado y todas las otras columnas en otro. SI es la columna que se separa antes, a una altura de \(\pm\) 25 en el dendograma, donde se separa de SPI y SGI, l que refuerza la disimilaridad de SI con todas las otras columnas, algo observado en ejercicios previos.
  • La separación que sucede a continuación de la ya mencionada es la que separa NS, CL y CT de PLI, CH y CI.
  • A su vez, NS se separa de CT y CL, mientras que PLI se separa de CI y CH.
  • Finalmente, solo quedan clústeres con columnas individuales o con grupos de dos columnas. CT y CL se separan a una altura de 14, SPI y SGI a una altura de 10 y CH y CI a una altura de 2, siendo las columnas más cercanas, algo que ya habíamos visto con las matrices de distancias.

Más o menos, estudiando el dendograma manualmente (algo que no es eficaç cuando hay muchas columnas), podemos establecer un orden de cercanía entre columnas/conglomerados.

  1. (CI) y (CH).
  2. (SPI) y (SPG).
  3. (CL) y (CT).
  4. (CI+CH) y (PLI).
  5. (CL+CT) y (NS).
  6. (CI+CH+PLI) y (CL+CT+NS).
  7. (SPI+SPG) y (SI).
  8. (CI+CH+PLI+CL+CT+NS) y (SPI+SPG+SI).

Además, para visualizar un análisis de conglomerados jerárquico también podemos usar un mapa de calor o heatmap. Para esto, podemos usar la función heatmap() del paquete stats, así como heatmap.2() (gplots), pheatmap() (pheatmap) o Heatmap() (ComplexHeatmap)

gplots::heatmap.2(hclust_1g$merge %>% as.matrix())
<b>Figure 1.7.</b>

Figure 1.7.

H

El siguiente paso es estudiar por algún criterio el número óptimo de conglomerados para el análisis jerárquico. Con la distancia de Canberra según R en particular, lo más sencillo es utilizar el criterio de las siluetas.

Las siluetas son un método para comprobar la consistencia de los conglomerados, de manera que representan como cada observación se parece a su cluster y al cluster más cercano. El valor de la silueta va de -1 a 1, con 1 siendo el valor que indica que más se parece a su propio cluster y -1 el que menos.

Cuando comparamos distintos numeros de clusters, podemos usar el valor medio de la silueta de cada numero de conglomerados para así elegir el número óptimo.

avsi <- c(0)
for(i in 2:8){
  hclust_1g <- hcut(x = dist_canb, k = i, isdiss = T, hc_method = "ward.D2")
  avsi[i-1] <- hclust_1g$silinfo$avg.width
}
plot(2:8,avsi,type="o", xlab="Number of clusters", ylab="Average silhouette length", main = "Hierarchical clustering")
<b>Figure 1.8.</b>

Figure 1.8.

Como se observa en la figura 1.8, hay dos picos en 4 clústeres (pico local) y en 6 clústers (pico absoluto). Con esto podemos decir que el número òptimo de conglomerados es 6.

I)

Estudiar con la misma distancia el número óptimo de conglomerados con el método PAM.

avsi <- c(0)
for(i in 2:8){
  km_vars <- pam(dist_canb, k=i, diss = T)
  si <- silhouette(km_vars, dist_canb)
  avsi[i-1] <- mean(si[,"sil_width"])
}
plot(2:8,avsi,type="o",xlab="Number of clusters", ylab="Average silhouette length", main = "PAM (K-medioids)")
<b>Figure 1.9.</b>

Figure 1.9.

Curiosamente, usando el método PAM( (figura 1.9), podmeos ver como el pico de la amplitud media de la silueta en los 6 conglomerados desaparace, convirtiéndose en un “valle”, mientras que el pico en los 4 clústeres sigue estando y el pico absoluto está en los 7 grupos.

En este caso, el número óptimo sería los de 7 grupos, aunqué en ambos, los 4 conglomerados parece un número razonable.

avsi <- c(0)
for(i in 2:8){
  km_vars <- pam(dist_canb, k=i, diss = T)
  si <- silhouette(km_vars, dist_canb)
  avsi[i-1] <- mean(si[,"sil_width"])
}
plot(2:8,avsi,type="o",xlab="Number of clusters", ylab="Average silhouette length", main = "PAM (K-medioids)")
<b>Figure 1.10.</b>

Figure 1.10.

J)

De los apartados anteriores se deduce que hay un número razonable de conglomerados, aunque no sea óptimo. Dibujar el dendograma del apartado (g) con esa partición.

hclust_1j <- hcut(x = dist_canb, k = 4, hc_method = "ward.D2", graph = T)
fviz_dend(hclust_1j, cex = 2 , k = 4, main = "Hierarchical clustering with Canberra distances", sub = "4 clusters",
          ggtheme = theme_bw(base_size = 17), show_labels = T, rect = TRUE)
<b>Figure 1.11.</b>

Figure 1.11.

Ejercicio 2.

Un estudio contiene dos medidas de los anillos de crecimiento en la escala del salmón de Alaska y de Canadá. Los datos se pueden obtener en el libro de Johnson et al. y se adjuntan en el archivo salmon.txt.

salmon <- read.delim("salmon.txt", sep = " ") %>% mutate(Gender = factor(Gender))
salmon_canada <- salmon %>% filter(Origin == "Canadian")
salmon_alaska <- salmon %>% filter(Origin == "Alaskan")

salmon_list <- list(salmon, salmon_canada, salmon_alaska) %>% set_names(nm = c("All salmons", "Canadian salmons", "Alaskan salmons"))
salmon_list_num <- salmon_list %>% purrr::map(~dplyr::select(.data = .x, Freshwater, Marine))

A)

Realizar una estadística descriptiva univariante y multivariante según el factor Origin. Añadir algunos gráficos ilustrativos, en particular el de dispersión.

  • Resumen descriptivo de los datos (extremos, quartiles, media…): El resumen descriptivo de los datos mediante summary() nos ayuda a obtener información sobre como están estructurados los valores de los datos (centralidad, etc). En este caso vamos a hacerlo en la totalidad de los datos, así como separados por origen.
purrr::map(salmon_list, summary)
## $`All salmons`
##  Gender   Freshwater        Marine           Origin  
##  1:52   Min.   : 53.0   Min.   :301.0   Alaskan :50  
##  2:48   1st Qu.: 99.0   1st Qu.:367.0   Canadian:50  
##         Median :117.5   Median :396.5                
##         Mean   :117.9   Mean   :398.1                
##         3rd Qu.:140.0   3rd Qu.:428.2                
##         Max.   :179.0   Max.   :511.0                
## 
## $`Canadian salmons`
##  Gender   Freshwater        Marine           Origin  
##  1:26   Min.   : 90.0   Min.   :301.0   Alaskan : 0  
##  2:24   1st Qu.:124.0   1st Qu.:345.2   Canadian:50  
##         Median :140.0   Median :369.5                
##         Mean   :137.5   Mean   :366.6                
##         3rd Qu.:151.5   3rd Qu.:385.8                
##         Max.   :179.0   Max.   :438.0                
## 
## $`Alaskan salmons`
##  Gender   Freshwater         Marine           Origin  
##  1:26   Min.   : 53.00   Min.   :355.0   Alaskan :50  
##  2:24   1st Qu.: 86.25   1st Qu.:402.0   Canadian: 0  
##         Median : 99.00   Median :427.5                
##         Mean   : 98.38   Mean   :429.7                
##         3rd Qu.:109.00   3rd Qu.:452.5                
##         Max.   :131.00   Max.   :511.0
  • Estadísticos descriptivos univariantes (media, mediana, desviación estándard) según origen: Con la media y la mediana tenemos información sobre la centralidad de los datos, mientras que con la desviación estándard obtenemos datos sobre la variabilidad de estos respecto la media. En este caso, vamos a hacerlo respecto a cada origen (salmón canadiense vs. salmón de Alaska.)
salmon %>% group_by(Origin) %>% summarise(`Mean | Freshwater` = mean(Freshwater),
                                          `Mean | Marine` = mean(Marine),
                                          `SD | Freshwater` = sd(Freshwater),
                                          `SD | Marine` = sd(Marine),
                                          `Median | Freshwater` = median(Freshwater),
                                          `Median | Marine` = median(Marine))
## # A tibble: 2 x 7
##   Origin `Mean | Freshwa… `Mean | Marine` `SD | Freshwate… `SD | Marine`
##   <fct>             <dbl>           <dbl>            <dbl>         <dbl>
## 1 Alask…             98.4            430.             16.1          37.4
## 2 Canad…            137.             367.             18.1          29.9
## # … with 2 more variables: `Median | Freshwater` <dbl>, `Median | Marine` <dbl>

Algo que es fácil de ver es que los salmones Alaskan suelen ser mayores que los Canadian cuando son marinos, pero menores en tamaño cuando se trata de salmones de agua dulce.

  • Matriz de varianzas-covarianzas:
var_alaska <- salmon %>% filter(Origin == "Alaskan") %>% dplyr::select(-Gender, -Origin) %>% var()
var_canada <- salmon %>% filter(Origin == "Canadian") %>% dplyr::select(-Gender, -Origin) %>% var()

var_alaska %>% knitr::kable(caption = "Variance-Covariance matrix for Alaskan salmons")
var_canada %>% knitr::kable(caption = "Variance-Covariance matrix for Canadian salmons")
Variance-Covariance matrix for Alaskan salmons
Freshwater Marine
Freshwater 260.6078 -188.0927
Marine -188.0927 1399.0861
Variance-Covariance matrix for Canadian salmons
Freshwater Marine
Freshwater 326.0902 133.5049
Marine 133.5049 893.2608
  • Varianza total:
paste("Varianza total de los salmones de Alaska ", sum(eigen(var_alaska)$values) %>% round(2))
paste("Varianza total de los salmones de Canada ", sum(eigen(var_canada)$values) %>% round(2))
## [1] "Varianza total de los salmones de Alaska  1659.69"
## [1] "Varianza total de los salmones de Canada  1219.35"
  • Varianza generalizada:
paste("Varianza generalizada de los salmones de Alaska ", det(var_alaska) %>% round(2))
paste("Varianza generalizada de los salmones de Canada ", det(var_canada) %>% round(2))
## [1] "Varianza generalizada de los salmones de Alaska  329233.85"
## [1] "Varianza generalizada de los salmones de Canada  273460.04"
  • Coeficiente de dependencia global:
paste("Coeficiente de dependencia global  de los salmones de Alaska ", 
      1-det(cor(salmon %>% filter(Origin == "Alaskan") %>% dplyr::select(-Gender, -Origin) )) %>% round(2))
paste("Coeficiente de dependencia global de los salmones de Canada ", 
      1-det(cor(salmon %>% filter(Origin == "Canadian") %>% dplyr::select(-Gender, -Origin) )) %>% round(2))
## [1] "Coeficiente de dependencia global  de los salmones de Alaska  0.1"
## [1] "Coeficiente de dependencia global de los salmones de Canada  0.0600000000000001"
  • Gráfico de dispersión: nos permitiran observar la correlación entre variables (Freshwater vs Marine) así como entre grupos (Canadian vs Alaskan).
car::scatterplot(salmon_alaska$Freshwater, salmon_alaska$Marine, 
                   main = "Alaskan Freshwater Salmons vs Alaskan Marine Salmons", xlab="Freshwater", ylab = "Marine")
Figure 2.1. Alaskan Freshwater Salmons vs Alaskan Marine Salmons

Figure 2.1. Alaskan Freshwater Salmons vs Alaskan Marine Salmons

car::scatterplot(salmon_canada$Freshwater, salmon_canada$Marine, 
                   main = "Canadian Freshwater Salmons vs Canadian Marine Salmons", xlab="Freshwater", ylab = "Marine")
Figure 2.2. Canadian Freshwater Salmons vs Canadian Marine Salmons

Figure 2.2. Canadian Freshwater Salmons vs Canadian Marine Salmons

car::scatterplot(salmon_alaska$Freshwater, salmon_canada$Freshwater, 
                   main = "Alaskan Freshwater Salmons vs Canadian Freshwater Salmons", xlab="Alaskan", ylab = "Canadian")
Figure 2.3. Alaskan Freshwater Salmons vs Canadian Freshwater Salmons

Figure 2.3. Alaskan Freshwater Salmons vs Canadian Freshwater Salmons

car::scatterplot(salmon_alaska$Marine, salmon_canada$Marine, 
                   main = "Alaskan Freshwater Marine vs Canadian Marine Salmons", xlab="Alaskan", ylab = "Canadian")
Figure 2.4. Alaskan Freshwater Marine vs Canadian Marine Salmons

Figure 2.4. Alaskan Freshwater Marine vs Canadian Marine Salmons

Vemos como ni los salmones del mismo lugar (ex. Alaska o Canada) tienen correlación entre agua dulce y agua marina (figuras 2.1 y 2.2), así como los del mismo tipo de agua tienen correlación entre los distintos lugares de origen (figuras 2.3 y 2.4).

  • Gráfico de cajas
salmon %>% melt() %>% 
  ggplot(aes(Origin, value, color = Origin)) + geom_boxplot() + facet_wrap(~variable, scales = "free") + 
  theme_pubr(legend = "none") +
  
  ylab("Growth") + xlab("Origin")
Figure 2.4. Boxplot of growth of freshwater and Marine salmons from Alaska and Canada.

Figure 2.4. Boxplot of growth of freshwater and Marine salmons from Alaska and Canada.

B)

Realizar un análisis discriminante lineal. La función lda() del paquete MASS puede servir.

Un análisis discriminante es un método usado para predecir la probabilidad de pertenecer a un grupo determinado basándonos en una o varias variables predictoras. En parte, es similar a la regresión logística (vista en otro tema) en el sentido que ambos métodos son usados para classificar datos en base a variables discretas (ej. grupo 1 y grupo 2), aunqué la regresión logística se usa en casos de variables discretas de dos clases, mientras que el análisis discriminante puede servir para 2 o más clases. 4 También podríamos usar más de una variable como variable predictora.

Cabe destacar que se necesita un grupo de datos conocidos para “entrenar” el modelo y así poder predecir la clasificación de datos nuevos.

Primero hay que usar la función lda() para calcular el modelo, luego, mediante la función predict(), para predecir las probabilidades de cada observación en pertenecer en un grupo u otro (tabla 2.1), así como los coeficientes del model lineal y la predicción de la clase, normalmente según el criterio de probabilidad > 0.5.

En este caso, la variable que defina los grupos será Origin, aunqué podríamos hacerla con otra variable categorica (ej. Gender) o crear una nueva variable que mezcle ambas.

salmon2 <- salmon %>% mutate(group = paste(Origin, Gender, sep = "_"))
lda_2b <- lda(salmon[,2:3], grouping = salmon$Origin)

lda_2b_pred <- lda_2b %>% predict()
lda_2b_pred$posterior %>% round(2) %>% head %>% knitr::kable(caption = "Table 2.1. Probabilities of belonging to Alaskan or Canadian categories for each of the salmons of the dataset. Note that only the first 6 observations are taken.")
Table 2.1. Probabilities of belonging to Alaskan or Canadian categories for each of the salmons of the dataset. Note that only the first 6 observations are taken.
Alaskan Canadian
0.43 0.57
0.02 0.98
1.00 0.00
1.00 0.00
0.93 0.07
0.99 0.01

A cada predicción se le asigna una puntuación o score relativo al modelo lineal que se relaciona con la probabilidad de pertenecer a uno u otro grupo y que nos puede dar información sobre la calidad de las predicciones.

ggplot(salmon, aes(lda_2b_pred$x, fill = Origin)) + geom_density(alpha = 0.5) + xlab("scores") 
<b>Figure 2.12. Density of the prediction scores in each group (Alaska / Canada) </b>

Figure 2.12. Density of the prediction scores in each group (Alaska / Canada)

Como se observa en la figura 2.12, hay un punto del gráfico de densidades en las que los valores de los dos grupos se solapan (cerca del 0), lo que nos puede sugerir que esas predicciones pueden no ser fiables, aunque en su mayor parte, los scores parece buenos.

Otra manera de mirar la calidad del modelo es usando una tabla de clasifiación cruzada, mostrando cuantos de cuandos son predecidos como un grupo son acertados y cuantos no.

table(lda=lda_2b_pred$class,real=salmon$Origin) %>%
  knitr::kable(caption = "Table 2.2. Table of cross-classification showing the predictied classes and the real ones.")
Table 2.2. Table of cross-classification showing the predictied classes and the real ones.
Alaskan Canadian
Alaskan 44 1
Canadian 6 49

Vemos en la tabla 2.2 que, de los 50 salmones canadienses, 49 han sido predichos como tal, mientras que en el caso de los salmones de Alaska han sido predichos correctamente en 44 casos, mostrando un poco menos de tasa de acierto, pero igualmente con un resultado aceptable.

C)

Clasificar una observación con un Freshwater de 120 y un valor de Marine de 380.

lda_2b <- lda(formula = Origin ~ Freshwater + Marine, data = salmon)

predict(lda_2b,newdata=data.frame(Freshwater=120,Marine=380))
## $class
## [1] Canadian
## Levels: Alaskan Canadian
## 
## $posterior
##     Alaskan  Canadian
## 1 0.2298261 0.7701739
## 
## $x
##         LD1
## 1 0.4199577

La predicción es que el salmón es canadiense, mientras que los salmones de medidas parecidas en los datos originales son tanto de Alaska como de Canadá.

D)

Comparar las matrices de covarianzas de las dos poblaciones con el test de la razón de verosimilitudes. También se puede aplicar el test M de Box. Ambos son muy sensibles a la no normalidad de los datos y tienden a rechazar la igualdad de covarianzas.

El análisis discriminante tiene asunciones, una de las cuales es la igualdad entre las mastrices de covarianzas entre las poblaciones que se están comparando. Para comprobar esto, podemos usar un test de verosimilitudes o un test M de Box.

Como sabemos que ambos test son sensibles a la no normalidad de los datos, primero haremos un test de shapiro multivariante para comprobar si estos son normales

mshapiro_test(data = salmon[,2:3])
## # A tibble: 1 x 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.995   0.970

Vemos como el p-valor del test de Shapiro-Wilk multivariante es muy proximo a uno, por lo que no podemos rechazar la hipótesis nula de una distribución normal, por lo que podemos fiarnos del resultado de los test.

  • Test de verosimilitudes:

  • Test M de Box: el test M de Box, o prueba de box permite comparar las matrices de covarianzas entre dos grupos para ver si estas son iguales. Es un test muy sensible y debido a esto, se recomienda un límite de significancia bastante bajo (i.e. 0.001). También hay que tener en cuenta, que en caso de falta de normalidad de los datos, el resultado del test de Box puede ser significativo debido a esta falta de normalidad y no a la diferencia entre las matrices de covarianzas [explicación sacada de la respuesta de mi PEC1]. En nuestro caso, la distribución parece ser normal, por lo que no debe preocuparnos esta parte. Para llevar a cabo el test, se puede usar la función boxM() del paquete heplots.

boxM_res <- heplots::boxM(salmon[,2:3], salmon$Origin, )
summary(boxM_res)
## Summary for Box's M-test of Equality of Covariance Matrices
## 
## Chi-Sq:   10.69615 
## df:   3 
## p-value: 0.01349 
## 
## log of Covariance determinants:
##  Alaskan Canadian   pooled 
## 12.70452 12.51891 12.72333 
## 
## Eigenvalues:
##     Alaskan Canadian    pooled
## 1 1429.3568 923.1148 1147.0461
## 2  230.3371 296.2362  292.4764
## 
## Statistics based on eigenvalues:
##               Alaskan    Canadian      pooled
## product   329233.8474 273460.0441 335483.8619
## sum         1659.6939   1219.3510   1439.5224
## precision    198.3702    224.2669    233.0522
## max         1429.3568    923.1148   1147.0461

Como se puede observar, tanto con la prueba de razón de verosimilitudes que con la prueba de Box, obtenemos que las matrices de covarianzas de los dos grupos son significativamente distintas bajo un nivel de significancia del 95% (p valor inferior a 0.05) y como la muestra parece seguir una distribución normal, esta diferencia es fiable.

Sin embargo, el test M de Box es muy sensible, por lo que se recomienda un límite de significancia de 0.001, por lo que, como el p-valor no es inferior a este, no podemos rechazar con total seguridad la hipótesis nula de covarianzas iguales y, por tanto, el resultado del análisis discriminante lineal es fiable.

E)

En el caso de poblaciones normales con diferentes matrices de covarianzas se clasificará cada observación en el grupo con máxima probabilidad a posteriori, pero entonces las funciones discriminantes no son lineales, ya que tienen un término de segundo grado. Realizar un análisis discriminante cuadrático. La función qda() del paquete MASS nos ayudará

qda_2e <- qda(formula = Origin ~ Freshwater + Marine, data = salmon)
qda_2e_pred <- predict(lda_2b)
table(qda=qda_2e_pred$class,real=salmon$Origin) %>%
  knitr::kable(caption = "Table 2.3. Table of cross-classification showing the predictied classes using the quadratic discriminant analysis and the real ones.")
Table 2.3. Table of cross-classification showing the predictied classes using the quadratic discriminant analysis and the real ones.
Alaskan Canadian
Alaskan 44 1
Canadian 6 49

En la tabla 2.3 se observan las predicciones de las observaciones originales y se comparan con el grupo original de dichas observaciones. El resultado es el mismo que en la tabla 2.2 pese a que se trata de un análisis discriminante lineal, mientras que en este caso estamos observando el resultado de un análsisi discriminante cuadrático.

F)

Calcular el número de parámetros que hay que estimar en la discriminación lineal y en la cuadrática.

El número de parámetros efectivos a estimar en el LDA se puede calcular con la siguiente función:

\[ K_p + (K-1) \]

Donde \(K\) es el numero de clases de la variable predictora.

G)

Calcular los errores de clasificación con ambas reglas utilizando validación cruzada. Si son similares, nos quedaremos con el análisis lineal que además es más robusto y de mejor interpretación.

Si vemos la tabla 2.2, podemos ver la tabla de datos predichos originales del análisis lineal, mientras que la tabla 2.3 muestra la misma para el análisis cuadrático.

table_lda <- table(lda=lda_2b_pred$class,real=salmon$Origin)
table_lda %>% knitr::kable(caption = "Table 2.2 (rep)")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_lda))/sum(table_lda),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_lda))/sum(table_lda)),2),"%")
Table 2.2 (rep)
Alaskan Canadian
Alaskan 44 1
Canadian 6 49
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"
table_qda <- table(qda=qda_2e_pred$class,real=salmon$Origin)
table_qda %>% knitr::kable(caption = "Table 2.3. (rep)")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_qda))/sum(table_qda),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_qda))/sum(table_qda)),2),"%")
Table 2.3. (rep)
Alaskan Canadian
Alaskan 44 1
Canadian 6 49
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"

Vemos como ambas tablas son iguales.

La validación cruzada es un método que nos permite evaluar el resultado de un análisis estadístico y que consiste en repetir el análisis mediante la evaluación de diferentes particiones (grupos de los datos originales que nos permiten entrenar el modelo).

Para llevar a cabo la validación cruzada podemos hacerlo mediante las funciones lda() y qda() indicando el argumento CV = T.

  • Validación cruzada para el LDA
lda_2b_CV <- lda(salmon[,2:3], grouping = salmon$Origin, CV = T)

table_lda_cv <- table(prediction = lda_2b_CV$class, real=salmon$Origin)
table_lda_cv %>% knitr::kable(caption = "Table 2.4. Table of predicted vs original data using LDA with CV.")

paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_lda_cv))/sum(table_lda_cv),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_lda_cv))/sum(table_lda_cv)),2),"%")
Table 2.4. Table of predicted vs original data using LDA with CV.
Alaskan Canadian
Alaskan 44 1
Canadian 6 49
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"
  • Validación cruzada para el QDA
qda_2h_CV <- qda(salmon[,2:3], grouping = salmon$Origin, CV = T)

table_qda_cv <- table(prediction = qda_2h_CV$class, real=salmon$Origin)
table_qda_cv %>% knitr::kable(caption = "Table 2.5. Table of predicted vs original data using QDA with CV.")

paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_qda_cv))/sum(table_qda_cv),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_qda_cv))/sum(table_qda_cv)),2),"%")
Table 2.5. Table of predicted vs original data using QDA with CV.
Alaskan Canadian
Alaskan 45 3
Canadian 5 47
## [1] "Percentage of coincidences without cross-validation: 92%"
## [1] "Error rate without cross-validation:: 8%"

Vemos como, después de usar la validación curzada (CV), los valores de acierto y error de los datos predecidos son muy similares a antes de usar la validación cruzada. Sin embargo, antes de la CV ambos métodos dieron el mismo resultado en cuanto a tasa de acierto y error en las predicciones, pero después pese a ser muy similares, el método LDA parece tener una mejor actuación, por lo que nos quedamos con este.

Bibliografía

  1. Alboukadel Kassambara (2017). Correspondence Analysis: Theory and Practice. STHDA. Online aquí.

  2. Michael Greenacre (2008). La pràctica del anàlisis de correspondencias. Capítulo 4: Distancia ji-cuadrado e inercia. Fundación BBVA. Online aquí.

  3. Alboukadel Kassambara (2017). CA in R Using FactoMineR: Quick Scripts and Videos. STHDA. Online aquí.

  4. Alboukadel Kassambara (2018). Discriminant Analysis Essentials in R. STHDA. Online aquí.